% To do likelihood ratio test
clear all;
clc;
disp('To do likelihood ratio test_MK2');

%load MK2 data (Both tests based on same data)
load ('MK2_TP06YRJfinR5G2501000kMR7.mat','PathCountTotal','PathCountAvgProb','NUMBER_OF_YEAR');
% load PathCountTotal_MK2_TP06YRJfinR5G1000kMR7.mat;
% load PathCountAvgProb_MK2_TP06YRJfinR5G1000kMR7.mat;
filename='TP06_Y_RJ_fin_R5_G250_1000k_MR7_Likelihood Ratio Test_MK1_MK2.xls';

% take the avg data from PathCountAvgProb (MK2).
N2(:,:,:) = PathCountTotal(1,:,:,:);
P2(:,:,:) = PathCountAvgProb(1,:,:,:);
ttlData = sum(N2(:));


% to calculate log (LMK1) hyphothesis
N1(:,:) = sum(N2(1:7,:,:));
P1=zeros(7,8);
for i=1:7
    for j=1:8
        P1(i,j)=N1(i,j)./sum(N1(i,1:8));
    end
end


LLMK1=0;
for i=1:7
    for j=1:8     
            if N1(i,j)~=0
                Y =log(P1(i,j)).* N1(i,j)./NUMBER_OF_YEAR;
                LLMK1=LLMK1+Y;
            else
                LLMK1=LLMK1+0;
            end 
       
    end
end


% to calculate log (LMK2) hyphothesis
LLMK2=0;
for i=1:7
    for j=1:7
        for k=1:8

            if N2(i,j,k)~=0
  
                Y =log(P2(i,j,k)).* N2(i,j,k)./NUMBER_OF_YEAR;
                LLMK2=LLMK2+Y;
            else
                LLMK2=LLMK2+0;
            end 
        end
    end
end


%Likelihood Ratio Test
invalid_Num_N1 = prod(size(find(N1==0)));
invalid_Num_N2 = prod(size(find(N2==0)));
dof = (prod(size(P2))-invalid_Num_N2) - (prod(size(P1))-invalid_Num_N1);% minus the number of zero element in P1 and P2.

[LRh,LRp,LRstat,cV] = lratiotest(LLMK2,LLMK1,dof);


% PP2(:,:)=P2(2,:,:);
% NN2(:,:)=N2(2,:,:); 
% sheetlist1=strcat('Prob of ij');
% sheetlist2=strcat('Prob of ijk_pathR2');
% sheetlist3=strcat('ObservationNum N1');
% sheetlist4=strcat('ObservationNum N2_pathR2');
% sheetlist5=strcat('Parameters');
% % xlswrite(filename, P1, sheetlist1);
% % xlswrite(filename, PP2, sheetlist2);
% % xlswrite(filename, N1, sheetlist3);
% % xlswrite(filename, NN2, sheetlist4);
% d = {'LRh','LRp','LRstat','cV','dof','LLMK1','LLMK2','ttlData'; LRh, LRp, LRstat,cV,dof,LLMK1,LLMK2,ttlData};
xlswrite(filename, d, sheetlist5)


